#include "fluid.h"

#include <omp.h>

#include <GL\glew.h>

#define GLFW_NO_GLU
#include <gl\glfw.h>

#include <cmath>
#include <basecode\misc_utils.h>
#include "basecode\shader.h"

float roygbiv[256][3];
float num = 12;

//#define SET_GRID( x, y, z, p ) grid[ x + (y*grid_width) + (z*grid_height*grid_width) ] = p; 
//#define GET_GRID( x, y, z ) grid[ x + (y*grid_width) + (z*grid_height*grid_width) ]; 
// x + (y*xSize) + (z*ySize*xSize)

	FluidSolver::FluidSolver()
	{
		N = (int) (num*num*num); //1024 * 2;
		G_factor = 0.02f;					// Gravity vec is unit length, scale with this constant
		G = glm::vec3(0.f,0.0,0.0f);		// gravity vector
		spacing = 3.0;						// Spacing of particles
		k = (real) (spacing / 1000.0);		// Far pressure weight
		k_near = k*10;						// Near pressure weight
		rest_density = 15;					// Rest Density
		sigma = 0.004f;						// visco 1
		beta = 0.004f;						// visco 2
		r = (real) (spacing*10.0);			// Radius of Support
		rsq = r*r;							// ... squared for performance stuff
		inv_r = 1.0f / r;                   // used to divide up grid and calc particle grid index
		SIM_W         = 300;				// The size of the world
		bottom        = 0;					// The floor of the world

		attracting = 0;
		repelling  = 0;
		attractor.x = -SIM_W;
		attractor.y = -SIM_W;
		do_viscosity_step = false;

		grid_width = (int)ceil( SIM_W * inv_r );
		grid_height= (int)ceil( SIM_W * inv_r );
		grid_depth = (int)ceil( SIM_W * inv_r );
		grid = new particle*[ grid_width * grid_height* grid_depth ]; 
		
		for(int i=0; i<grid_width * grid_height* grid_depth; i++)
		{
			grid[i] = NULL;
		}


		// x + (y*xSize) + (z*ySize*xSize)
		initParticles();
		generateColorLUT();

		glGenBuffers(1, &VertexVBOID);
		glBindBuffer(GL_ARRAY_BUFFER, VertexVBOID);
		glBufferData(GL_ARRAY_BUFFER, verts.size()*sizeof(glm::vec3), NULL, GL_DYNAMIC_DRAW); // GL_STREAM_DRAW - per render update, GL_DYNAMIC_DRAW per frame

		colorShader = new Shader();
		colorShader->addStage(std::pair<GLenum, const GLchar*>(GL_VERTEX_SHADER, "color.vert"));
		colorShader->addStage(std::pair<GLenum, const GLchar*>(GL_FRAGMENT_SHADER, "color.frag"));
		//colorShader->install();
		colorShader->end();
	}

	void DrawCube(float size)
	{
		// method taken from cube.c by Mark J. Kilgard, 1997
		// http://www.opengl.org/resources/code/samples/glut_examples/examples/cube.c
		float n[6][3] = {  // Normals for the 6 faces of a cube. 
			{-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 0.0, 0.0},
			{0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, 0.0, -1.0} };

			int faceIndices[6][4] = {  /* Vertex indices for the 6 faces of a cube. */
				{0, 1, 2, 3}, {3, 2, 6, 7}, {7, 6, 5, 4},
				{4, 5, 1, 0}, {5, 6, 2, 1}, {7, 4, 0, 3} };


				float v[8][3];  // Will be filled in with X,Y,Z vertexes.
				v[0][0] = v[1][0] = v[2][0] = v[3][0] = -size * 0.5f;
				v[4][0] = v[5][0] = v[6][0] = v[7][0] = size * 0.5f;
				v[0][1] = v[1][1] = v[4][1] = v[5][1] = -size * 0.5f;
				v[2][1] = v[3][1] = v[6][1] = v[7][1] = size * 0.5f;
				v[0][2] = v[3][2] = v[4][2] = v[7][2] = size * 0.5f;
				v[1][2] = v[2][2] = v[5][2] = v[6][2] = -size * 0.5f;

				glBegin(GL_QUADS);
				for (int i=0; i < 6; i++) {
					glNormal3fv(&n[i][0]);
					glVertex3fv(&v[faceIndices[i][0]][0]);
					glVertex3fv(&v[faceIndices[i][1]][0]);
					glVertex3fv(&v[faceIndices[i][2]][0]);
					glVertex3fv(&v[faceIndices[i][3]][0]);
				}
				glEnd();

	}


	void FluidSolver::render()
	{
		static float _minVelocity = 0;
		static float _maxVelocity = 0;
		float mindens = -rest_density;
		
#define USE_POINT_SPRITE
#ifdef USE_POINT_SPRITE
		glEnable(GL_TEXTURE_2D);
		glEnable(GL_BLEND);
		glEnable(GL_POINT_SPRITE);
		glPointSize(r*2);
		
		glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
		//glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_CONSTANT_COLOR); // used with glBlendColorEXT, asserts default state is glBlendColor(0,0,0,1);
		glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
		glDepthMask(GL_FALSE);
#else
		glDisable(GL_BLEND);
		glDisable(GL_TEXTURE_2D);
		glDisable(GL_POINT_SPRITE);
		glPointSize(1);
#endif


		if( particlesCreated > 0 )
		{
			double before_upload = glfwGetTime();
//#define USE_VBO
#ifdef USE_VBO
			verts.clear();
			// TODO, might get better perf if positions are seperate from particles... and just upload positions directly to vbo...meh.
			for(int i=0; i<particlesCreated; ++i){
				verts.push_back( particles[i].pos );
			}

			//glEnableClientState(GL_VERTEX_ARRAY);
			//glVertexPointer( 3, GL_FLOAT, 0, &verts[0] );
			//glDrawArrays(GL_POINTS, 0, verts.size() );
			//glDisableClientState(GL_VERTEX_ARRAY);

			glBindBuffer( GL_ARRAY_BUFFER, VertexVBOID );
			glEnableClientState(GL_VERTEX_ARRAY);
			#define BUFFER_OFFSET(i) ((char *)NULL + (i))
			glVertexPointer(3, GL_FLOAT, sizeof(glm::vec3), BUFFER_OFFSET(0));   //The starting point of the VBO, for the vertices

			//int attrib_loc = glGetAttribLocation(colorShader->getProgram(), "vVtxPosition");
			//glVertexAttribPointer(attrib_loc, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), &verts[0].x );
			//glEnableVertexAttribArray(weightPosition);

			glBufferData(GL_ARRAY_BUFFER, verts.size()*sizeof(glm::vec3), &verts[0].x, GL_DYNAMIC_DRAW); // Use insead of sub, cause num particles changes...
			//glBufferSubData(GL_ARRAY_BUFFER, 0, verts.size()*sizeof(glm::vec3), &verts[0].x);
			glDrawArrays(GL_POINTS, 0, verts.size() );
			glDisableClientState(GL_VERTEX_ARRAY);
#else
					glBegin(GL_POINTS);
					for(int i=0; i < particlesCreated; ++i)
					{
			//#define velocitycolor
			#ifdef velocitycolor
						float velMagnitude = glm::length( particles[i].vel );
						if( velMagnitude < _minVelocity )
						{
							_minVelocity = velMagnitude;
						}
						if( velMagnitude > _maxVelocity )
						{
							_maxVelocity = velMagnitude * 0.86f;
						}
						float finalC = ( (_minVelocity+velMagnitude)/_maxVelocity  );
						int c = std::min( std::max( (int) (170.0f - finalC*170.0f), 0 ), 170 );
						glColor3fv( roygbiv[c] ); // same as glColor3f( roygbiv[c][0], roygbiv[c][1], roygbiv[c][2] );	
			#else  // Use pressure for color

						float opacity = ( particles[i].rho ); // try log10()
						//float f = 2.0f * particles[i].force.len();
						
						opacity/=(rest_density-mindens);
						float finalC = opacity*3.0f;
						
						//float finalC = opacity*opacity + (3 - 2*opacity);
						//glColor4f( f, f, opacity, opacity);
						int c = std::min( std::max( (int) (170.0f - finalC*170.0f), 0 ), 170 );
						//glColor3fv( roygbiv[c] ); // same as 
						glColor4f( roygbiv[c][0], roygbiv[c][1], roygbiv[c][2], opacity );		
			#endif
						glVertex3fv( &particles[i].pos.x );
					}
					// Done drawing points
					glEnd();
#endif


			uploadTime = glfwGetTime() - before_upload;

#pragma region draw_cube
			glEnable(GL_DEPTH_TEST);
			glDisable(GL_POINT_SPRITE);
			glDisable(GL_BLEND);
			glDisable(GL_TEXTURE_2D);
			glDepthMask(GL_TRUE);
			glColor3f(1,1,1);

			glPushMatrix();
			glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
			glTranslatef( SIM_W*0.5f, SIM_W*0.5f, SIM_W*0.5f );
			DrawCube( SIM_W );
			glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
			glPopMatrix();
#pragma endregion draw_cube
			//glPushMatrix();
			//glTranslatef( SIM_W*0.5f, SIM_W*0.5f, SIM_W*0.5f );
			//glm::vec3 grav_arrow = G;
			//glm::normalize(grav_arrow);
			//grav_arrow *=  150.0f;
			//glColor3f(0,1,0);
			//glBegin(GL_LINES);
			//glVertex3f(0,0,0);
			//glVertex3fv( &grav_arrow.x );
			//glEnd();
			//glPopMatrix();
		}
	}


	void FluidSolver::update()
	{
		advance_particles_verlet();
		update_grid();
		calc_density();
		calc_pressure();
		if (do_viscosity_step) calc_viscosity();
}

void FluidSolver::advance_particles_euler(){
	// UPDATE
	//
	// This modified verlet integrator has dt = 1 and calculates the velocity
	// For later use in the simulation. 

	// For each particles i ...	
#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i)
	{			
		// Apply the currently accumulated forces
		particles[i].vel += particles[i].force;
		particles[i].pos += particles[i].vel;
		particles[i].vel *= 0.9995f;

		// Restart the forces with gravity only. We'll add the rest later.
		particles[i].force = G * G_factor;

		// If the particle is outside the bounds of the world, then
		// Make a little spring force to push it back in.
		real rdist = r * 0.5f;
		if(particles[i].pos.x < rdist) particles[i].force.x -= (particles[i].pos.x-rdist ) / 2;
		if(particles[i].pos.x >  SIM_W-rdist) particles[i].force.x -= (particles[i].pos.x - (SIM_W-rdist)) / 2;
		if(particles[i].pos.y < rdist ) particles[i].force.y -= (particles[i].pos.y - rdist ) / 2;
		if(particles[i].pos.y > SIM_W-rdist) particles[i].force.y -= (particles[i].pos.y - SIM_W+rdist) / 2;
		if(particles[i].pos.z < rdist) particles[i].force.z -= (particles[i].pos.z-rdist ) / 2;
		if(particles[i].pos.z >  SIM_W-rdist) particles[i].force.z -= (particles[i].pos.z - (SIM_W-rdist)) / 2;

		// Handle mouse interaction. 

		glm::vec3 p_to_attr =  (particles[i].pos - attractor);
		real attr_dist2 = glm::dot(p_to_attr, p_to_attr); // Dot with self is the same as dist squared
		real attr_l = 32;

		//  make mouse attract particles
		if( attracting )
			if( attr_dist2 < attr_l*attr_l && attr_dist2 > 16 ){
				glm::normalize(p_to_attr);
				particles[i].force -= p_to_attr * 0.1f;// * (1.0f / 64.0f );
			}


			//  make mouse repel particles
			if( repelling )
				if( attr_dist2 < attr_l*attr_l ){
					particles[i].force += (particles[i].pos - attractor) * (1.0f / 64.0f);
				}

				// Reset the nessecary items.
				particles[i].rho = 0;
				particles[i].rho_near = 0;
	}
}

void FluidSolver::advance_particles_verlet(){
	// UPDATE
	//
	// This modified verlet integrator has dt = 1 and calculates the velocity
	// For later use in the simulation. 

	// For each particles i ...	
	#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i)
	{			
		// Apply the currently accumulated forces
		particles[i].pos += particles[i].force;

		// Restart the forces with gravity only. We'll add the rest later.
		particles[i].force = G * G_factor;

		// Calculate the velocity for later.
		particles[i].vel = particles[i].pos - particles[i].pos_old;

		// Normal verlet stuff
		particles[i].pos_old = particles[i].pos;
		particles[i].pos += particles[i].vel;

		// If the particle is outside the bounds of the world, then
		// Make a little spring force to push it back in.
		real rdist = r * 0.5f;
		if(particles[i].pos.x < rdist) particles[i].force.x -= (particles[i].pos.x-rdist ) / 2;
		if(particles[i].pos.x >  SIM_W-rdist) particles[i].force.x -= (particles[i].pos.x - (SIM_W-rdist)) / 2;
		if(particles[i].pos.y < rdist ) particles[i].force.y -= (particles[i].pos.y - rdist ) / 2;
		if(particles[i].pos.y > SIM_W-rdist) particles[i].force.y -= (particles[i].pos.y - SIM_W+rdist) / 2;
		if(particles[i].pos.z < rdist) particles[i].force.z -= (particles[i].pos.z-rdist ) / 2;
		if(particles[i].pos.z >  SIM_W-rdist) particles[i].force.z -= (particles[i].pos.z - (SIM_W-rdist)) / 2;

		// Handle mouse interaction. 

		glm::vec3 p_to_attr =  (particles[i].pos - attractor);
		real attr_dist2 = glm::dot(p_to_attr, p_to_attr); // Dot with self is the same as dist squared
		real attr_l = 32;

		//  make mouse attract particles
		if( attracting )
			if( attr_dist2 < attr_l*attr_l && attr_dist2 > 16 ){
				glm::normalize(p_to_attr);
				particles[i].force -= p_to_attr * 0.1f;// * (1.0f / 64.0f );
			}


			//  make mouse repel particles
			if( repelling )
				if( attr_dist2 < attr_l*attr_l ){
					particles[i].force += (particles[i].pos - attractor) * (1.0f / 64.0f);
				}

				// Reset the nessecary items.
				particles[i].rho = 0;
				particles[i].rho_near = 0;
	}
}


void FluidSolver::update_grid()
{
	#pragma omp parallel for
	for(int x=0; x<grid_width; x++){
		for(int y=0; y<grid_height; y++){
			for(int z=0; z<grid_depth; z++){
				//SET_GRID(x,y,z,NULL);
				grid[ x + (y*grid_width) + (z*grid_height*grid_width) ] = NULL; 
			}
		}
	}
	#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i){
		particle *p = &particles[i];
		int gridx = (int) floor( p->pos.x * inv_r );
		int gridy = (int) floor( p->pos.y * inv_r );
		int gridz = (int) floor( p->pos.z * inv_r );
		particle *successor = grid[ gridx + (gridy*grid_width) + (gridz*grid_height*grid_width) ]; //GET_GRID(gridx,gridy,gridz);
		//grid[gridx][gridy][gridz] = p;
		//SET_GRID(gridx,gridy,gridz,p);
		grid[ gridx + (gridy*grid_width) + (gridz*grid_height*grid_width) ] = p; 
		p->succ = successor;
	}
}

void FluidSolver::calc_density()
{
	// Calculate the density by basically making a weighted sum
	// of the distances of neighboring particles within the radius of support (r)
	const int dirx[] = { -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1 };
	const int diry[] = { -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1 };
	const int dirz[] = { -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
	#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);
		int int_z = (int) floor(particles[i].pos.z * inv_r);

		for(int j=0; j<27; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];
			int gridz = int_z + dirz[j];

			if( gridx < 0 || gridy < 0 || gridz < 0 ||  gridx >= grid_width || gridy >= grid_height || gridz >= grid_depth ) continue;

			particle *p2 = grid[ gridx + (gridy*grid_width) + (gridz*grid_height*grid_width) ];

			while( p2 != NULL )
			{
#pragma region CALC_DENSITY
				if( p1->id > p2->id ) // dont check self or same twice
				{
					// The vector seperating the two particles
					glm::vec3 rij = p2->pos - p1->pos;

					// Along with the squared distance between
					real rij_magnitudesquared = glm::dot( rij, rij );         

					// If they're within the radius of support ...
					if(rij_magnitudesquared < rsq)
					{
						// Get the actual distance from the squared distance.
						real rij_len = sqrtf(rij_magnitudesquared);

						// And calculated the weighted distance values
						real q = (real) 1.0f - rij_len / r;
						real q2 = q*q;
						real q3 = q2*q;

						p2->rho += q2;
						p2->rho_near += q3;
						particles[i].rho      += q2;
						particles[i].rho_near += q3;

					} // eof in range
				} // eof is >
#pragma endregion CALC_DENSITY
				p2 = p2->succ;
			} // while!null
		} // for grid

	} // for each particle...
}


void FluidSolver::calc_pressure()
{
	// PRESSURE
	// Make the simple pressure calculation from the equation of state.
	#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i)
	{
		particles[i].press = k * (particles[i].rho - rest_density);
		particles[i].press_near = k_near * particles[i].rho_near;
	}

	// PRESSURE FORCE
	// We will force particles in or out from their neighbors
	// based on their difference from the rest density.
	const int dirx[] = { -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1 };
	const int diry[] = { -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1 };
	const int dirz[] = { -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
	#pragma omp parallel for
	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);
		int int_z = (int) floor(particles[i].pos.z * inv_r);

		for(int j=0; j<27; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];
			int gridz = int_z + dirz[j];

			if( gridx < 0 || gridy < 0 || gridz < 0 || gridx >= grid_width || gridy >= grid_height || gridz >= grid_depth ) continue;

			particle *p2 = grid[ gridx + (gridy*grid_width) + (gridz*grid_height*grid_width) ];
			while( p2 != NULL )
			{
				if( p1->id > p2->id ) // only compare particles once. also makes sure we dont compare same. of course.
				{
					glm::vec3 rij = p2->pos - p1->pos;
					real rij_magnitudesquared = glm::dot( rij, rij );         
					
					if(rij_magnitudesquared < rsq)
					{
						real rij_len = sqrtf(rij_magnitudesquared);

						// And calculated the weighted distance values
						real q = (real) 1.0 - rij_len / r;
						real q2 = q*q;

						// calculate the force from the pressures calculated above
						real dm = (p1->press + p2->press) * q +
							(p1->press_near + p2->press_near) * q2;

						// Get the direction of the force
						//vec3f D = rij.normal() * dm;

						real forceMag = (real) 1.0/rij_len * dm;
						glm::vec3 D = rij * forceMag;    //  OPTIMISATION!! here we save a sqrt() call.
						//dX += D;
						p2->force += D;
						p1->force -= D;
					} // eof in range
				}
				p2 = p2->succ;
			}
		} // eof grid loop
		//p1->force -= dX;
	} // eof force-loop
}
void FluidSolver::calc_viscosity()
{
	const int dirx[] = { -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1 };
	const int diry[] = { -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1 };
	const int dirz[] = { -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 };


	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);
		int int_z = (int) floor(particles[i].pos.z * inv_r);

		for(int j=0; j<27; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];
			int gridz = int_z + dirz[j];

			if( gridx < 0 || gridy < 0 || gridz < 0 || gridx >= grid_width || gridy >= grid_height || gridz >= grid_depth ) continue;

			particle *p2 = grid[ gridx + (gridy*grid_width) + (gridz*grid_height*grid_width) ];
			while( p2 != NULL )
			{         
				if( p1->id > p2->id ) // only compare particles once. also makes sure we dont compare same. of course.
				{
					glm::vec3 rij = p2->pos - p1->pos;
					real len = sqrt( glm::dot(rij,rij) );
					real q = len / r;
	     
					glm::vec3 rij_normal = rij / len;
					 // Get the projection of the velocities onto the vector between them.
					glm::vec3 vel = p2->vel - p1->vel; 
					real u = glm::dot( rij_normal, vel );
	         
					 if(u > 0.0f)
					 {
						// Calculate the viscosity impulse between the two particles
						// based on the quadratic function of projected length.
						real magnitude = (1.0f - q) * (p2->sigma * u + p2->beta * u*u);
						glm::vec3 I = rij_normal * magnitude;
			       
						// Apply the impulses on the two particles
						/*
						particles[ni_idx].vel.x -= I_x * 0.5;
						particles[ni_idx].vel.y -= I_y * 0.5;
			            
						particles[nj_idx].vel.x += I_x * 0.5;
						particles[nj_idx].vel.y += I_y * 0.5;
						*/
			            
						// Using verlet, so move positions forward instead with force.
						p1->force -= (I * 0.5f);
						p2->force += (I * 0.5f);
					 }
				}

				p2 = p2->succ;
			} // while( p2 != NULL )
		}
	}
}



void FluidSolver::activateAttractor( int mx, int my )
{
	attracting = 1;
	attractor = glm::vec3( (real )mx, -my + SIM_W, 0.0f);	
}

void FluidSolver::activateRepelling( int mx, int my )
{
	repelling = 1;
	attractor = glm::vec3( (real )mx, -my + SIM_W, 0.0f);	
}

void FluidSolver::deactivateAttractor()
{
	attracting = 0;
	repelling = 0;
}

void FluidSolver::addParticles( int mx, int my )
{
	attractor.x = (real) SIM_W*0.5f;
	attractor.y = (real) SIM_W*0.75f;	
	attractor.z = (real) SIM_W*0.5f;	
	if (particlesCreated < MAX_PARTICLES)
	{
		real x = attractor.x + MISRENDER::frand() * 20.0f;
		real y = attractor.y + MISRENDER::frand() * 20.0f;
		real z = attractor.z + MISRENDER::frand() * 20.0f;
		particles[particlesCreated].pos = glm::vec3(x,y,z);
		particles[particlesCreated].pos_old = glm::vec3(x,y,z);  
		particles[particlesCreated].force = glm::vec3(0.0f,0.0f,0.0f);
		particles[particlesCreated].sigma = sigma;
		particles[particlesCreated].beta = beta;
		particles[particlesCreated].id = particlesCreated;

		particlesCreated++;
	}
}

void FluidSolver::initParticles()
{
	// create a world with dimensions x:[-SIM_W,SIM_W] and y:[0,SIM_W*2]
	// Initialize particles
	// We will make a block of particles with a total width of 1/2 of the screen.
	particlesCreated = 0;

	real startvalue = SIM_W * 0.25f;
	real x = startvalue;
	real y = startvalue;
	real z = startvalue;
	real limit_pos = SIM_W * 0.75f; // Don't allow particles outside the box
	for(int i=0; i<N; i++)
	{
		if(particlesCreated > N) break;

		particles[particlesCreated].pos = glm::vec3(x,y,z);
		particles[particlesCreated].pos_old = glm::vec3(x,y,z);
		particles[particlesCreated].vel = glm::vec3(0.0f,0.0f,0.0f);
		particles[particlesCreated].force = glm::vec3(0.0f,0.0f,0.0f);
		particles[particlesCreated].sigma = sigma;
		particles[particlesCreated].beta = beta;
		particles[particlesCreated].id = particlesCreated;

		particlesCreated++;		

		real inc = (SIM_W*0.5f) / num;
		x += inc;
		if( x >= limit_pos)
		{ 
			x = startvalue;
			y += inc;
		}

		if( y >= limit_pos)
		{
			y = startvalue;
			z += inc;
		}

		if( z >= limit_pos)
		{
			z = SIM_W * 0.75f;
		}

	}
	verts.resize( N );
}

void FluidSolver::generateColorLUT()
{
	// Generate color look-up table
	double hue = 0;
	double hueInc = 6.0/255.0;
	double v = 1.0; // value 
	double s = 1.0; // saturation
	for(int i=0; i<175; i++)
	{
		double c = v * s;
		double x = c * 1.0-fabs(fmod(hue,2.0)-1.0);
		//x*=255;
		//c*=255;
#define RGB(i, r,g,b) roygbiv[i][0] = (float)r; roygbiv[i][1] = (float)g; roygbiv[i][2] = (float)b;
		if(hue >= 0 && hue <= 1){ RGB(i, c,x,0); }// red->orange
		if(hue >= 1 && hue <= 2){ RGB(i, x,c,0); }// orange->yellow
		if(hue >= 2 && hue <= 3){ RGB(i, 0,c,x); }// yellow->green
		if(hue >= 3 && hue <= 4){ RGB(i, 0,x,c); }// green->blue
		if(hue >= 4 && hue <= 5){ RGB(i, x,0,c); }// blue->indigo
		if(hue >= 5 && hue <= 6){ RGB(i, c,0,x); }// indigo->violet
#undef RGB
		hue += hueInc;
	}
}

FluidSolver::~FluidSolver()
{
	glDeleteBuffers(1, &VertexVBOID );
	delete[] grid;
}

void FluidSolver::reinstallShader()
{
	colorShader->install();
}
